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I.  INTRODUCTION 


Filamentary  organic  matrix  composites,  hereafter  called  composites,  are 
efficient  structural  material,  coupling  low  weight  with  high  strength.  They 
have  been  used  for  many  years  in  various  secondary  and  primary  structural 
roles  both  in  commercial  products  and-  in  military  applications.  A  primary  use 
is  for  pressure  containment  structures  such  as  rocket  motor  cases,  launch 
tubes,  and  cold  gas  storage  vessels.  Some  design  philosophies  require  that 
such  pressure  containment  structures  be  lOOZ  proof  tested,  i.e.,  proof  test 
each  item  produced.  These  proof  tests,  to  be  effective  in  eliminating  defec¬ 
tive  hardware  and  assuring  safety  during  the  operational  life,  are  usually 
conducted  at  significant  load  levels.  It  is  a  well-knotra  and  frequently  ob¬ 
served  phenomenon  that  such  proof  tests  induce  a  form  of  damage  called  crazing 
in  most  composites.  Other  forms  of  damage  may  also  occur  such  as  matrix/fiber 
debonding  and  ply  delamination.  Hereafter  the  term  ‘‘damage”  will  be  used  to 
denote  the  internal  damage  resulting  from  a  proof-test  loading.  In  order  to 
fully  utilize  the  high  specific  strength  of  composites,  a  mechanism  to  include 
the  effect  of  this  proof-test-induced  damage  is  needed. 

There  has  been  a  significant  amount  of  work  in  the  area  of  composite  dam¬ 
age  with  a  large  part  being  devoted  to  fatigue  damage.  A  compilation  of  some 
of  these  efforts  can  be  found  in  References  [1],  [2],  [3],  [4],  and  [5].  The 
following  discussion  %d.ll  point  out  some  of  the  different  approaches  used  to 
include  damage  accumulation  in  materials. 

Smith  and  Huang  [6]  and  Lee  [7]  have  used  linearly  elastic  finite  element 
analyses  to  model  damage  accumulation  in  composites.  Both  use  a  failure  cri¬ 
terion  to  define  the  matrix  damage  zone.  The  stiffness  is  then  reduced  in  the 
damaged  elements.  The  loading  is  incremental  until  a  global  failure  condition 
is  met. 

Relfsnider  and  Hlghsmlth  [8],  [9]  have  investigated  stiffness  reduction  in 
laminated  composites.  Tliey  reference  an  earlier  work  by  Relfsnider  [10]  where 
a  characteristic  damage  state  of  crack  patterns  in  transverse  plies  was  pre¬ 
dicted.  The  crack  patterns  were  shown  to  be  a  laminate  property,  determined 
only  by  lamina  properties  and  the  orientations  and  stacking  of  the  laminae. 

In  References  [8]  and  [9],  a  linear  fracture  mechanics  approach  was  taken  to 
predict  crack  growth  and  subsequent  failure  where  the  moduli  were  reduced  as  a 
function  of  the  crack  geometry. 

Nulsmer  and  Tan  [II]  included  the  effects  of  matrix  cracking  through  a 
constitutive  theory  based  upon  simple  physical  models  of  ply  damage,  where 
certain  parameters  relating  to  the  growth  of  ply  damage  and  subsequent  loss  of 
ply  stiffness  are  predicted  using  a  Griffith-type  energy  balance.  Several 
Investigators  [12-15]  have  used  what  they  termed  "continuum  damage  mechanics” 
to  study  the  accumulation  of  damage  in  an  elastic  homogeneous  body.  This 
approach  is  based  on  earlier  work  by  Kachanov  [16]  where  the  process  of  fail¬ 
ure  was  considered  as  a  process  of  crack  formation  against  a  background  of 
creep  deformation.  Kachanov  defines  a  quantity  called  the  "continuity”  which 
is  a  measure  of  the  damage  or  cracking  in  the  body.  Murakami  [12]  elucidates 
the  notion  of  continuum  damage  mechanics  which  hypothesizes  that  the  effects 
of  damage  to  a  homogeneous  body  can  be  described  by  appropriate  mechanical 
variables  »rtilch  he  termed  damage  variables.  The  continuity  function  defined 
by  Kachanov  is  extended  to  tensor  form.  Controlled  experiments  are  conducted 
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iflth  good  agreement  between  results  and  theory.  Krajclnovlc  [13,  15]  and 
Krajclnovlc  and  Fonseka  [14]  have  also  Investigated  damage  In  homogeneous 
elastic  materials  using  continuum  damage  mechanics  with  similar  success.  All 
of  the  efforts  discussed  used  either  a  linear  elastic  or  linear  fracture  me¬ 
chanics  approach  to  Incorporate  the  Internal  damage  of  a  material.  Schapery 
[17]  approached  the  problem  of  damage  In  composites  from  a  viscoelastic  point 
of  view.  Schapery  developed  relationships  of  the  hereditary  Integral  form  to 
Include  not  only  damage,  but  also  a  “healing, ^  l.e.,  decrease  In  damage, 
effect. 

In  view  of  the  well-known  viscoelastic  behavior  of  organic  matrix  compos¬ 
ites,  It  seems  appropriate  that  this  behavior  be  Included  In  the  simulation 
model  discussed  here.  Well-established  classical  viscoelastic  procedures  are 
used  with  modifications  to  account  for  damage  and  Its  accumulation.  The 
Inclusion  of  damage  will  be  shown  to  produce  a  nonlinear  model. 

As  noted  earlier,  Kachanov  [16]  used  a  parameter  termed  the  continuity 
ftmctlon  to  define  the  state  of  partial  cracking,  l.e.,  partial  damage,  of  a 
material.  Appealing  to  this  Idea,  a  damage  function,  u.  Is  defined.  It  Is 
assumed  that  this  function  has  zero  value  when  no  damage  Is  present  and 
Increases  monotonlcally  to  the  value  of  one  at  failure.  One  may  think  of  u  as 
an  Index  of  the  Internal  damage. 

A  viscoelastic  material  may  Initially  contain  minute  voids  or  other  de¬ 
fects.  As  the  material  Is  loaded  and  creep  occurs,  these  voids  will  tend  to 
grow  in  size.  New  voids  will  form  at  the  defects  and  grow  also.  At  some 
time,  these  voids  will  coalesce  Into  microcracks  and  continue  to  grow  as, a 
result  of  the  creep  of  the  material  until  critical  cracks  are  formed  and 
failure  occurs.  Other  forms  of  damage  may  also  be  aggravated  by  the  creep. 

The  creep  of  this  damaged  material  produces  a  change  In  Its  deformation  state. 
It  Is  well  known  that  a  change  In  displacement  state  Is  related  to  a  change  in 
the  state  of  strain  In  the  material.  This  Implies  that  strain  and  the  damage 
that  has  accumulated  re  related.  Therefore,  It  Is  assumed  that  the  damage  Is 
some  function  of  the  state  of  strain  that  exists. 

O'Brien  [18]  and  Hlghsmlth  and  Relfsnlder  [9]  have  shown  that  damage  accu¬ 
mulation  reduces  the  stiffness  of  a  composite.  Hlghsmlth  and  Relfsnlder 
modeled  the  overall  laminate  response  by  reducing  the  transverse  modulus  of 
the  laminae  as  a  result  of  accumulated  damage.  This  Idea  Is  central  to  the 
development  here  and  Is  the  basis  for  the  model  that  Is  constructed.  Here,  as 
In  Hlghsmlth's  and  Relfsnlder' s  work,  a  laminate  analysis  using  the  damage- 
reduced  transverse  modulus  of  the  laminae  Is  employed  to  predict  the  laminate 
stiffness  values.  These  resulting  stiffnesses  are  then  used  In  a  simulation 
model  to  predict  the  viscoelastic  response  of  a  composite  cylinder  which 
Includes  the  damage  effect.  The  simulation  uses  both  linear  elastic  and 
linear  viscoelastic  responses  modified  to  account  for  the  dependence  of  the 
transverse  modulus  on  the  strain  via  the  damage  function.  The  matrix  is 
assumed  to  be  linearly  viscoelastic  and  Isotropic  with  the  fiber  being 
linearly  elastic  and  Isotropic.  The  well-known  elastic-viscoelastic 
correspondence  principle  is  employed  in  the  solution.  A  discussion  of  the 
correspondence  principle  is  found  in  Section  V.  Basically,  It  allows  one  to 
obtain  the  Laplace  transform  of  a  time-dependent  solution  by  simple  substitu¬ 
tions  of  various  Laplace  transformed  parameters  into  an  elastic  solution  of 
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the  same  geometry  and  constraints.  This  Laplace  transformed  solution  Is 
called  the  "associated  elastic"  solution.  Due  to  the  complexity  of  the  asso¬ 
ciated  elastic  solution,  a  closed-form  Inversion  Is  not  sought.  Rather,  the 
associated  elastic  solution  Is  Inverted  by  a  numerical  procedure  to  yield  the 
time-dependent  response.  The  numerical  procedure  requires  that  the  associated 
elastic  solution  be  known  only  for  discreet  values  of  the  Laplace  parameter. 

The  remaining  sections  are  arranged  In  the  order  that  the  respective  ele¬ 
ments  are  used  In  the  solution  process.  Section  II  develops  the  governing 
equations  for  an  elastic  solution  for  a  cylindrical  body  under  uniform  Inter¬ 
nal  pressure.  These  equations  are  then  solved  to  yield  the  radial  displace¬ 
ment,  meridional  strain,  radial  strain,  and  circumferential  strain.  Section 
III  develops  the  elastic  material  relationships  property  that  Is  required  In 
the  solution.  The  transformation  from  local  lamina  to  global  laminate  proper¬ 
ties  Is  defined.  The  various  micromechanical  models  for  the  lamina  constitu¬ 
tive  properties  are  also  Included  for  completeness.  In  Section  IV,  a  short 
discussion  of  the  time-dependent  material  properties  Is  Included.  Section  V 
discusses  the  viscoelastic  solution  procedure.  The  correspondence  principle 
Is  defined  and  the  parameters  for  substitution  Into  the  elastic  solution  to 
produce  the  associated  elastic  solution  are  noted.  The  numerical  Inversion' 
technique  Is  discussed.  The  approximate  solutions  for  both  the  linearly 
varying  load  and  a  step  load  are  developed  and  the  method  for  determining  the 
unknown  solution  constants  are  shown.  Section  VI  contains  a  discussion  of 
damage  In  a  viscoelastic  material.  In  Section  VII,  Implementation  of  the 
solution  Is  discussed.  The  results  of  a  numerical  example  are  Included  In 
Section  VIII  and  a  listing  of  the  program  developed  Is  contained  In  the 
Appendix. 
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where  the  e^j  are  tensor  strain  components  and  the  are  displacements  In  the 
B,r,z  coordinate  system*  For  uniform  loading,  the  shear  strain  components 
vanish,  and  *  0.  At  a  point  sufficiently  far  from  the  ends  of  the 
cylinder,  e^z  assumed  to  be  constant.  Thus  Bqs.  (3)  reduce  to: 
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C.  Constitutive  Equations 


For  a  uniformly  loaded  cylinder,  the  shear  stresses  and  shear 
strains  are  zero  and  the  constitutive  relationship  reduces  to  [19]: 
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The  elements  of  the  3x3  matrix  of  Eq.  (5)  will  be  developed  in  a 
following  section. 

D.  Boundary  Conditions 

Boundary  conditions  for  the  internally  pressurized  cylinder  are: 

Ogj.  ■  -P  (at  the  inner  surface)  ( 

Opj,  •  0  (at  the  outer  surface)  ( 

where  P  is  the  internal  pressure. 

E.  Solution 

Combining  Eqs.  (2),  (4)  and  (5)  yields  the  governing  differential 

equation: 


The  subscript,  r,  has  been  dropped  on  the  radial  displacement  for  conven¬ 
ience,  l.e.,  ■■  u. 


The  general  solution  to  Eq.  (7)  Is  given  as  [20]: 
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The  parameter  a  Is  a  measure  of  the  anisotropy  In  the  plane  of  cross  section. 
Normally,  a  ■  1  only  for  unidirectional  composites  with  fiber  direction 
parallel  to  the  meridional  axis.  For  this  model,  the  fibers  are  not  parallel 
to  the  meridional  axis.  Therefore,  the  expression  given  by  Eq.  (9a)  will  be 
the  only  one  considered. 

From  Eqs.  (4b)  and  (9a): 
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From  Eqs.  (4a)  and  (9a): 
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Combining  Eqs.  (6),  (10)  and  (11)  yields: 
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Substituting  Eqs.  (6)  Into  Eq.  (16)  gives: 
-P  -  A(Cf^  +  oC„)ri"^^"“^  +  B(C^^  - 

■*■  *is  l^zr  ■*■  ^  ■*■ 


0  -  A(Cft.  +  +  BCCft.  -  oC„)ro‘^^'^“^ 

■*■  ®zz  ll'zr  ■•■  ^  ■‘'  (1®) 

Equations  (16),  (17)  and  (18)  can  now  be  solved  simultaneously  to  obtain  A,  B 
and  e^2 


-  (bii  -  b2i) 


-  -j  (Bi2  “  B22) 


A  ^'’13  “  '’23) 


where 


a22«33  ■  a32a23 

-(a2ia33  -  a3ia23) 
a2ia32  -  a22a3i 

-Cai2a33  ”  a32ai3) 

aiia33  -  a3iai3 

■(aiia32  *  a3iai2) 

aiia22a33  +  ai2a23a3i  +  ai3a2ia32 


-ai3a22a3i  -  ai2a2ia33  -  aiia32a23 


•  O  ^  a  •••  O  a*'  *  •'V  »»«-•, 


1  +  a 


■J  ^^2  fl 


(21a) 


«11 


?  L‘ 


*12 


h  [==^1  <■•- 


a  C,r) 


1 

2  2 


*13 

®  ^  0  +  * 

r2  **  1  - 

1 

*21  " 

<Cer  +  * 

*22  " 

(Cfir  -  «  C„)ri-(l+«) 

fl 

*23  " 

Czr  j  ^ 

*31  “ 

(Cqc  +  a  Crr)^o  ^ 

*32  - 

(C0r  -  Cl  Crr)ro“^^‘^®^ 

fl 

*33  - 

Czr  j  ^ 

(Cjfl  +  Cji-) 


(21b) 

(21c) 

(21d) 

(21e) 

(21f) 

(21g) 

(21h) 

(211) 


The  radial  dlaplacemenC  given  by  Eq.  (9a)  may  now  be  determined.  The  hoop 
strain,  cqs,  and  radial  strain,  crrt  >■■^7  be  calculated  by  using  Eqs.  (10)  and 
(11),  respectively.  With  the  strains  (crrt  £00  Ezz)  known,  the  stresses 
may  be  determined  from  Eq.  (5). 
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III.  MATERIAL  PROPERTIES 


The  cylinder  under  discussion  Is  asswsed  co  be  coaposed  of  a  lasdnate 
made  up  of  several  laminae.  It  Is  further  assumed  that  the  laminae  alternate 
with  an  angle  of  +  y  fo  the  hoop  direction  as  shown  In  Figure  1.  To  produce 
the  laminate  stiffness,  one  must  combine  the  stiffnesses  of  the  various  lami¬ 
nae  In  an  appropriate  manner.  The  first  step  Is  the  transformation  of  the 
lamina  properties  from  the  local  L,T,r  coordinate  system  to  the  global  9,s,r 
coordinate  system. 

A.  Lamina  Stiffness 


Consider  the  vector  transformation  from  the  lamina  L,T,r  coordinate 
system  to  the  global  9,z,r  system  given  by: 


(“el 


%  “1  ®i 


/  ^ 


«L 


< 


> 


^2  <>2  02 


(22 


I3  m3  03 


/ 


where  the  e^  are  the  respective  components  and  the  prime  denotes  the  global 
Q,t,z  system.  The  direction  cosines,  1^,  m;^,  n^^  (1  *  1,2,3)  are  given  by: 


• 

cos 

a,0) 

m 

cos  Y 

“1 

m 

cos 

(T,0) 

m 

sin  Y 

"1 

m 

cos 

(r,  9) 

m 

0 

m 

cos 

(L,z) 

m 

-sin  Y 

^2 

m 

cos 

(T,2) 

m 

cos  Y 

"2 

m 

cos 

(r,z) 

m 

0 

^3 

m 

cos 

<L,r) 

m 

0 

m 

cos 

(T,r) 

m 

0 

no 

m 

cos 

(r,r) 

m 

1 

(23a 

(23b 

(23c 

(23d 

(23e 

(23f 

(23g 

(23h 


(231 


It  l8  well  known  that  the  stresses  transform  according  to: 
rf  -  I  2 

where  ^  v 


2 

”l 

2 

"l 

2\mj 

2m^n3 

^”1*1 

2 

”2 

2 

"2 

CM 

B 

2m2n2 

^^2*2 

4 

2 

”3 

2 

"3 

2V3 

2111303 

203X3 

*l'2 

®1®2 

"l"2 

(  ^2  )  (®in2+njin2) 

*2*3 

"2"3 

(02  ^3+^' 

*3*1 

m3«^ 

"3"l 

(03X3+13 
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For  the  local  system,  the  constitutive  relationship  Is  given  by: 


2“  £  £  ( 
where  C  Is  a  6x6  matrix  of  effective  lamina  properties  and  Is  given  by  [21] 


^LL 

ClT 

CLr 

0 

0 

0 

^LT 

Ctt 

Ctr 

0 

0 

0 

CLr 

Cxr 

Crr 

0 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

0 

2  Gpi. 

0 

0 

0 

0 

0 

0 

2  G^r 

and  tensor  strain,  e,  given  by. 


A  single  bar  under  a  variable  Is  used  to  denote  a  vector  or  first-order 
tensor.  A  double  bar  under  a  variable  Is  used  to  denote  a  second-order  tensor 
quantity. 

The  effective  lamina  properties  are  given  by  [21]: 


ClL  “ 

Ommbi  S 

(1  -  El) 

6 

•  \t  '*rL> 

^TT 

6 

( MjfL  ^  '*rT^ 

^LT 

6 

Cfr  " 

^  '*rT  'tT'*rL^ 

6 

(31a) 


(31b) 


(31c) 


(31d) 


(1  -  \»rr  ■  (3ie) 

Gtt  -  - 26 - 

and,  because  of  Isotropy  In  the  T-r  plane: 

^Lr  “  ^LT  (31f) 

Cyj.  *  C<p<j'  (31g) 

^Lr  ”  (31h) 

with 

6  “  (1  +  '<]-r)(l  “  “  2  (311) 


The  constituent  property  models  for  El,  Ef,  vj-p,  vl^,  vjl  ^T  wlH 
discussed  later. 

Substituting  Eq.  (28)  Into  Eq.  (24): 

o'  -  T  C  e  (32) 


The  tensor  strain  transforms  according  to: 


Solving  for  e  In  Eq.  (33): 
e'  -  !“’■  €' 


where  Che  -1  superscript  denotes  Che  matrix  Inverse.  Substituting  Eq.  (34) 
Into  Eq.  (32),  and  noting  that  for  Che  primed,  global  coordinate  system 


C’  e' 


yields: 


T  C  T 


Thus  It  Is  seen  that 


C  -  T  C  T 


where  Is  the  stiffness  of  Che  lamina  transformed  to  Che  global  6,r,z 
system.  The  task  now  is  to  perform  Che  Indicated  multiplication  of  Eq.  (37) 
to  produce  the  C  values  In  Eq.  (5).  After  much  algebra,  one  obtains  for  the 
reduced  lamina  stiffness  transformed  to  the  global  6,r,z,  system: 

Cjj  ■  sln^Y  +  cos^Y  +  2  cos^y  sln^Y  (Clt  +  2  (38« 

C2  0  ”  ^LT  (slii^Y  cos^y)  "b  sln^  Y  COS^Y  (^LL  ^TT  ~  ^  ^HiT^  (381 

Cjj.  ■  sln^Y  Clj  +  cos^Y  Cxr  (38c 

CflQ  -  cos^Y  '^LL  Ctt  +  2  sin2Y  cos^Y  (Clt  +  2  Glt)  (38<i 

C  gr  ■  cos^Y  +  sln^Y  Cfr  (38« 


B.  Laminate  Stiffness 

Consider  now  a  laminate  ccmposed  of  several  lamina.  For  this  model, 
a  symmetric  laminate  with  N  laminae  of  alternating  +  y  wrap  angle  Is  assumed. 
Gamma  (  y)  is  shown  In  Figure  1  and  la  the  angle  between  the  fiber  direction, 
L,  and  the  global  circumferential  direction,  6.  For  the  kill  lamina  of  the 
laminate,  the  constitutive  relationship  Is  given  by; 


where  N  Is  the  total  number  of  laminae. 

Using  the  approach  taken  by  Chamls  [22],  the  reduced  properties  of  the 
laminate  may  be  found  from: 


1  ** 

■kE 


*^i  £i 


where  Is  the  3x3  matrix  In  Eq.  (39),  t^  Is  the  total  laminate  thickness, 
and  t^  Is  the  thickness  of  the  l^h  lamina.  From  Eqs.  (38),  It  Is  noted  that 
the  elements  of  are  even  functions  of  y  Equation  (40)  then  reduces  to: 


^zz  ^z8  ^zr 


^zfi  ^09  *^0r 
^zr  ^0r  ^rr 


C.  Constituent  Properties 

The  laminate  stiffness  has  now  been  expressed  In  terms  of  the  wrap 
angle  and  the  lamina  properties.  The  lamina  properties  In  turn  can  be 
expressed  In  terms  of  the  constituent  properties  and  volume  fractions  of  the 
composite. 

A  filamentary  composite  is  made  up  of  fibers  embedded  In  a  matrix. 
Both  the  fiber  and  the  matrix  are  assumed  herein  to  be  Isotropic.  The  longi¬ 
tudinal  modules,  E^,  and  the  major  Poisson's  ratio,  can  be  expressed  by  a 

rule  of  mixtures  formulation  [23]  as: 

El  -  VfEf  +  v^En,  (42) 

NLT  -  Vf  Vf  +  v^  (43) 

where  E^  Is  the  modulus  of  the  matrix,  Ef  Is  the  modulus  of  the  fiber,  Vf  is 
the  Poisson's  ratio  of  the  fiber,  is  the  Poisson's  ratio  of  the  matrix, 

Vf  Is  the  fiber  volume  fraction,  and  v^,  Is  the  matrix  volume  fraction. 

Expressions  for  vrj  are  Glx  af®  given  by  Whitney  [24]  and  Foye  [25] 
as,  respectively: 


vfvf  +  Vn,\te 


1  +  %  -  VLT  ^ 

1  -  4  +  %  vlT 


®m 

^LT  “  — 


4  -  It  +  nN 


(4-Tt)N  +  n 
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%ihere 


Gf(Ti  ^^4vf)  +  Gni(it  -4vf) 
Gf  (  It  -4vf  )  +  Gni(  It  +4vf  ) 


where  Gg^  Is  the  matrix  shear  modulus  and  Gf  is  the  fiber  shear  modulus. 


as: 


The  transverse  modulus,  Ef,  is  given  by  the  lower  bound  expression  [21] 


Ef 


1 

+  Ihl 

Ef  E„ 


(47 


From  a  consideration  of  the  83mimetry  of  the  compliance  matrix,  the  minor 
Poisson's  ratio,  vjxi  la  expressed  as: 


E'p 

'*rL  “  eT  'iT 


(48 


IV.  TIME-DEPENDENT  MATERIAL  PROPERTIES 


In  this  development.  It  Is  assumed  that  only  the  matrix  exhibits  time- 
dependent  properties,  these  being  designated  relaxation  modulus.  Eg,,  and 
Poisson's  ratio,  These  two  quantities  may  be  determined  from  experimental 

tests  of  the  matrix  material  In  three  different  ways  [26],  [27]. 

The  relaxation  modulus  may  be  measured  directly  In  a  relaxation  test  and 
a  method  such  as  least  squares  used  to  fit  the  data  to  a  Prony  series  of  the 
form: 


N 

E^(t)  -  Bj  exp(-Xjt)  +  E^ 

J-1  ' 


(49) 


where  N  Is  the  total  number  of  exponential  terms  Included  with  the  \j  being 
assigned.  The  Bj  and  E_  are  constants  determined  from  a  curve  fit  of  the 
data. 


A  second  method  to  determine  the  relaxation  modulus  Is  to  first  obtain, 
from  a  creep  test,  the  creep  compliance,  Jj,(t),  of  the  form: 

N 

J^(t)  -  G  +  Ht  +  ^  expC-R^t)  (50) 

1-1 


where  the  R^  are  assigned  and  with  the  P^,  and  H  being  constants  determined 
by  a  curve  fit  to  the  data.  An  expression  relating  the  Laplace  transforms  of 
the  creep  compliance  and  relaxation  modulus  [27]  Is: 


Em(8) 


1 

s^T  (s) 
tn 


(51) 


where  the  bar  over  a  variable  Is  used  to  denote  the  Laplace  transform  of  that 
function  and  s  In  the  Laplace  parameter.  Inversion  of  Eq.  (51)  yields  an 
expression  for.  the  relaxation  modulus  of  the  form  given  by  Eq.  (49).  A 
thorough  development  of  this  second  method  Is  given  by  Hackett  and  Dozier  In 
Reference  [26]. 

A  third  method  would  yield  the  relaxation  modulus  from  the  complex 
relaxation  modulus.  Details  of  this  method  are  given  by  Flndly  In  Reference 
[27]. 


The  time-dependent  Poisson's  ration  Is  obtained  from  the  well-known 
expression  for  the  bulk  modulus: 


Rearranging  Eq.  (52)  gives: 


Vt) 


1  Em(t) 

2  6K 


'7^^.-  -7-  ^  -7 


(53) 


where  the  bulk  modulus  Is  time-independent* 

V.  VISCOELASTIC  SOLUTION 

A.  Elastlc-Vlscoelastlc  Correspondence  Principle 

The  elastlc-vlscoelastlc  correspondence  principle,  hereafter  referred 
to  as  the  correspondence  principle.  Is  a  well-known  method  for  the  solution  of 
viscoelastic  problems  [26],  [29]  and  will  be  employed  here*  Basically,  the 
correspondence  principle  allows  one  to  convert  an  elastic  solution  for  a  given 
geometry  and  nonmoving  boundary  conditions  to  a  quasl-statlc  "associated 
elastic"  solution  In  the  Laplace  domain*  This  requires  that  the  corresponding 
substitutions  noted  In  Table  1  be  made  In  the  elastic  solution* 

TABLE  1*  Elastlc-Vlscoelastic  Correspondence 
Principle  Parameters 


Elastic 

Substitution 

Clj(t) 

j(8) 

qj(t) 

Hj(8) 

E„(t) 

8^m(®> 

VO 

8*^(8) 

Ti(t) 

Ti(s) 

All  elements  of  Table  1  have  been  defined  earlier  except  for  T{^(t)  which  Is 
the  time-varying  surface  traction*  This  associated  elastic  solution  Is  then 
inverted  back  Into  the  physical  domain  to  yield  the  time-dependent  solution. 

If  the  associated  elastic  solution  has  a  closed  form  expression.  It 
may  be  possible  to  use  an  exact  Inversion  technique,  e*g.,  the  method  of  par¬ 
tial  fractions,  to  obtain  the  time-dependent  solution.  However,  many  times 
the  solution  may  be  known  only  for  discrete  values  of  the  Laplace  parameter. 
If  this  Is  the  case,  then  some  numerical  Laplace  Inversion  technique  may  be 
used.  Even  in  cases  where  a  closed  form  solution  Is  possible,  the  complexity 
of  the  solution  may  dictate  the  use  of  a  numerical  Inversion  procedure*  Due 
to  the  complexity  of  this  problem,  the  latter  approach  Is  taken* 

B*  Numerical  Laplace  Transform  version 

There  are  several  numerical  Inversion  techniques  which  may  be  used 
[30]*  The  one  used  here  Is  called  the  method  of  collocation  [31]  and  Its 
application  has  been  demonstrated  by  numerous  researchers  [28],  [32],  [33]. 


*.K- 
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%■*  " 
« 


(54) 


Schapery  [31]  presented  a  general  expression  for  a  stress  or  displacement 
response,  4<C),  as: 

N 

(Kt)  -  4,’  +  (l,“t  (!iy  exp(-t/(i^) 

v-1 

where  (j;' ,  4;'*,  and  are  constants  to  be  determined,  are  prescribed  posi¬ 
tive  constants,  and  N  Is  the  number  of  exponential  terms  Included.  As  will  be 
discussed  later,  the  simulation  model  is  constructed  by  decomposing  the  pres¬ 
sure  load  Into  the  sum  of  several  loads  as  shown  In  Figure  2.  Superposition 
of  the  several  respective  responses  Is  then  used  to  determine  the  total  re¬ 
sponse.  From  Figure  2,  It  Is  noted  that  only  two  types  of  loadings  are  con¬ 
sidered:  a  linearly  varying  load  and  a  step  load. 

For  the  step  load  shown  In  Figure  3a,  It  has  been  shown  [31]  that  the 
circumferential  strain,  £90,  can  be  approximated  by: 

N 

€00(t)  -  A  hy  exp(-t/ay)  (55) 

v-1 

where  the  constant  A  la  given  by: 

N 

A  -  e'09(to)  -  ^  \  (56) 

v-1 

with  e00(tQ)  being  the  elastic  response  at  time  t^  and  with  the  prime  denoting 
a  step  load. 

Substituting  Eq.  (56)  Into  Eq.  (55),  taking  the  Laplace  transform  and 
rearranging  gives: 

8e09(to)  “  £00(8) 


where  eq@(s)  Is  the  associated  elastic  solution  for  the  circumferential  strain 
found  by  using  the  correspondence  principle  as  discussed  earlier  and  s  Is  the 
Laplace  parameter. 

For  this  model,  a  six-term  exponential  function  Is  used  to  define  the 
circumferential  strain  as  a  function  of  time,  l.e.,  N  -  6.  Therefore,  six 
discrete  values  of  the  Laplace  parameter,  s,  are  needed  to  solve  for  the  six 
unknown  constants,  h^.  The  choice  of  the  discrete  Laplace  parameters  Is  a 
matter  of  judgment  and  experience.  A  reasonable  expression  for  av  and  s  seems 
to  be: 


Off  -  exp(7-2v) 


(58a) 
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Utilizing  these  values,  Eq.  (57)  may  be  solved  for  the  h^'s.  The  constant  A 
is  then  calculated  using  Eq.  (56).  The  time-dependent  circumferential  strain 
may  then  be  found  from  Eq.  (55).  The  results  will  be  the  time-dependent  cir¬ 
cumferential  strain  due  to  a  stepwise  internal  pressure  load. 

The  second  type  of  load  to  be  considered  is  one  that  varies  linearly 
with  time.  It  is  assumed  that  the  circumferential  strain  may  be  expressed  by: 

N 

^15e(t)  -  C  +  Dt  gv  exp(-t/av)  (59) 

v-1 

where  the  double  prime  denotes  a  response  to  a  linearly  varying  load  with  C,  D 
and  gy  being  constants  to  be  determined. 

The  initial  conditions  for  the  linearly  varying  load  of  Figure  3b  are: 


*^efi(t) 

-  0 

at  t  -  0 

(6i)a) 

ae99(t) 

at 

eO0(ta) 

^a 

at  t  -  0 

(60b) 

where  e99(tg)  is 
arbitrary  time, 

the  elastic 
tg,  as  shown 

response  to  a 
in  Figure  3b. 

pressure,  P^,  corresponding 
The  time  derivative  of  Eq. 

to  some 
(59) 

is: 

N 

— -  -  —  exp(-t/0 

v-1 


Substituting  Eqs.  (60)  Into  Eq.  (59)  and  Eq*  (61)  and  rearranging  gives: 


Substituting  Eqs>  (62)  Into  Eq.  (59)  yields: 

"  4(!  (t.)  "  8, 

eesCt)  ■  -  8,  +  t+  ^  — 

V“1  V*1 


+  gv  exp(-t/av) 


Taking  the  Laplace  transform  of  Eq.  (63)  gives: 


E  «v  — S- 

s  + - 


where  eq0(s),  the  Laplace  transform  of  ep0(t).  Is  the  associated  elastic  solu¬ 
tion  determined  by  using  the  correspondence  principle  as  discussed  earlier.  A 
six-term  exponential  function  Is  used  to  define  the  circumferential  strain, 
l.e.,  N  >  6.  Thus,  six  arbitrary  values  of  the  Laplace  parameters,  s,  are 
needed  to  solve  for  the  six  unknowns,  g^.  This  choice  Is,  as  before,  given  by 
Eqs.  (58).  Substituting  Eq.  (58b)  Into  Eq.  (64)  and  rearranging  gives: 


/  °v«w  .  4  _  \ 

+  (V  Ov 


C9fi(s) 


E9fl(ta)  ^2 
— t -  ^ 


from  which  the  gy's  may  be  found.  The  constants,  C  and  D,  may  then  be 
obtained  using  Eqs.  (62).  Substitution  of  C  and  D  Into  Eq.  (59)  gives  the 
time-dependent  circumferential  strain  resulting  from  the  linearly  varying 
loads. 


v"*  .*•  »*■*  L’* 
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VI.  DAMAGE  IN  A  VISCOELASTIC  MATERIAL 


A  viscoelastic  material  may  Initially  contain  minute  voids  or  other 
defects.  As  the  material  Is  loaded  and  creep  occurs,  the  voids  will  grow  in 
size.  New  voids  will  form  at  the  defects  and  grow  also.  At  some  time,  these 
voids  will  coalesce  into  microcracks  and  continue  to  grow,  again  as  a  result 
of  the  creep  of  the  material.  This  condition  progresses  until  critical  cracks 
are  formed  and  failure  occurs.  Other  forms  of  damage  may  also  be  aggravated 
by  the  creep.  Kachanov  [11]  Introduced  a  parameter  termed  the  "continuity"  to 
account  for  the  damage  In  an  elastic  material.  Borrowing  this  Idea,  a  damage 
function,  iii,  Is  defined.  Assuming  the  damage  to  be  Isotropic,  the  damage 
function  becomes  a  scalar.  The  damage,  {4  can  be  thought  of  as  an  Index  of 
the  damage  that  exists  In  the  materials.  A  damage  value  of  zero  denotes  no 
damage.  The  value  1  is  arbitrarily  assigned  to  the  damage  at  failure.  The 
deformation  of  the  material  Is  related  to  the  creep  chat  occurs.  It  Is  well 
known  that  Che  strain  Is  a  function  of  the  deformation.  Therefore,  It  Is 
assumed  that  Che  damage  may  be  expressed  as  a  scalar  function  of  the  strain. 
Gerstle  [29],  using  the  maximum  strain  failure  theory,  has  demonstrated  good 
correlation  between  Che  circumferential  strain  and  failure  of  blaxlally-loaded 
quasl-lsotroplc  cylinders.  The  failure  strain,  Kq,  being  given  as,  for 
linearly  elastic  fibers: 


4Ef 


where  ?f  Is  the  fiber  fracture  strength.  In  view  of  this  correlation.  It  Is 
assumed,  as  a  first  approximation,  that  the  damage  Is  a  function,  of  the  cir¬ 
cumferential  strain  only.  Many  choices  exist  for  this  relationship.  Curve 
#1  In  Figure  4  represents  some  general  relationship  which  includes  a  "healing" 
or  decrease  of  damage  with  a  decrease  In  circumferential  strain.  A  first 
approximation  for  the  circumferential  strain-damage  relationship  Is  denoted  by 
Curve  #2  In  Figure  4,  where  the  damage  increases  monotonlcally  with  circum¬ 
ferential  strain.  This  results  In  Che  damage  remaining  constant  for  a 
decrease  In  circumferential  strain  as  shown  by  the  horizontal  portion  of  Curve 
#2  In  Figure  4.  The  corresponding  analytical  expression  for  the  linearly 
varying  portion  of  Curve  #2  Is: 
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where  e.  Is  the  circumferential  failure  strain  and  is  set  equal  to  K0  given  by 
Eq.  (68).  Failure  of  the  cylinder  will  then  occur  when  the  following 
expression  Is  satisfied: 

€gQ(t)  >  (68) 

It  can  be  enlightening  to  consider  the  response  of  a  simple  viscoelastic 
material.  A  schematic  response  diagram  for  a  simple  sprlng-dashpot  model  is 
shown  In  Figure  5.  The  load  shown  In  Figure  6  Is  applied  to  this  simple 
model.  A  sketch  of  Che  resultant  response  Is  shown  In  Figure  7.  The  straight 
lines  with  slopes  kj  and  kj  +  k2  are  the  respective  elastic  response  curves 
for  stiffnesses  kj  and  kj  +  k2  where  kj  and  k2  are  the  respective  spring 


CIRCUMFERENTIAL  STRAIN 


Figure  4.  Damage  versus  circumferential  strain 


constants  for  springs  1  and  2.  As  the  load  Increases  to  point  A,  the  response 
Is  defined  by  segment  OA  In  Figure  7.  At  point  A,  damage  occurs.  This  damage 
Is  simulated  by  cutting  one  of  the  springs,  say  k|^.  In  order  for  a  force 
balance  to  be  maintained,  the  force  In  the  dashpot  must  Increase.  This 
results  In  an  Increase  In  the  displacement  rate  which  Is  Indicated  by  the 
change  In  slope  of  the  response  curve  at  point  A.  The  load  continues  to 
Increase  until  point  B  Is  reached.  The  corresponding  response  Is  given  by 
segment  AB  In  Figure  7.  At  point  B,  the  load  Is  held  constant.  If  held 
constant  long  enough,  the  displacement  will  Increase  to  point  C  which  Is  an 
equilibrium  point  and  lies  on  the  elastic  response  curve,  k2.  However,  at 
point  C  the  load  Is  decreased  and  the  deformation  follows  some  curve  from 
point  C  to  point  D.  At  point  D,  the  load  Is  zero,  but  the  displacement  Is  not 
fully  recovered.  The  recovery  continues  along  segment  DO  with  decreasing 
speed  until  the  deformation  reaches  zero.  After  a  sufficiently  long  time,  the 
only  evidence  that  remains  of  the  damage  Is  the  reduced  stiffness.  The 
response  of  this  or  any  other  model  to  future  loading,  e.g..  In-service 
conditions,  must  be  based  upon  this  deduced  stiffness.  The  point  to  be  empha¬ 
sized  here  Is  that  damage  Is  modeled  as  a  reduction  In  stiffness. 

Several  experimental  studies  have  produced  data  which  relate  the  stiff¬ 
ness  of  a  composite  system  to  an  Increase  In  Internal  damage  [8],  [9],  [35], 
[36],  [37].  Hlghsmlth  and  Relfsnlder  [9]  modeled  the  overall  response  of  a 
laminate  containing  damage  by  reducing  the  transverse  modulus  of  the  laminae. 
O'Brien  [18]  used  the  change  In  transverse  modulus  as  a  predictor  of  buckling 
In  a  damaged  composite.  This  basic  Idea  of  reducing  the  transverse  modulus  of 
the  laminate  as  a  function  of  damage  Is  employed  here. 

Several  different  forms  for  the  transverse  modulus-damage  relationship 
are  possible.  Hlghsmlth  and  Relfsnlder  [9]  published  data  which  Indicates 
large  stiffness  changes  early  In  the  life  of  a  composite  material  undergoing 
fatigue  loading.  The  pertinent  data  Is  reproduced  here  In  Figure  8.  As  a 
first  approximation.  It  Is  assumed  that  a  transverse  modulus-damage  rela¬ 
tionship  Is  similar  to  the  stlffness/number-of-cycles  relationship  shown  In 
Figure  8.  A  simple  model  of  the  relationship  between  transverse  modulus  and 
damage  Is  shown  In  Figure  9.  The  corresponding  analytical  expression  Is: 

Et  -  (1  -  oj)2  for  0  <  1  (69) 

where  E-p  Is  the  Initial  or  undamaged  transverse  modulus.  Substituting  Eq. 

(67)  Into  Eq.  (69)  yields: 

Et  -  Et^  ^1  -  (70) 


where  the  transverse  modulus  now  becomes  a  function  of  time  as  a  consequence 
of  the  circumferential  strain  being  time-dependent.  In  addition,  the  for¬ 
mulated  problem  has  become  nonlinear  because  of  the  Interdependence  of  £99  and 
Et.  In  the  simulation  model,  this  nonlinearity  Is  accounted  for  by  an  Itera¬ 
tive  process  whereby,  the  first  transverse  modulus  Is  determined  from  Eq.  (70), 
using  the  circumferential  strain  from  the  previous  Iteration.  Convergence  Is 
based  upon  the  value  of  circumferential  strain  calculated  at  a  given  time  for 
the  Initial  linearly  varying  load  and  on  time  to  failure  for  the  remaining  por¬ 
tion  of  the  simulation. 


VII.  SIMULATION  MODEL 


In  this  section  the  strategy  used  to  develop  the  simulation  model  Is 
discussed.  The  simulation  model  developed  predicts  the  circumferential  strain 
of  a  filament-wound  cylinder  having  accumulated  Internal  damage.  The  proce¬ 
dure  followed  Is  to  first  decompose  the  load  Into  several  components  as  shown 
In  Figure  2.  The  respective  solutions  are  calculated  and  superposition  used 
to  construct  the  total  solution.  For  load  0  <  t  an  elastic  solution, 

which  Includes  the  damage  effect  given  by  Eq.  (69),  Is  computed.  This 
approach  Is  taken  for  this  phase  of  the  loading  since  there  Is  very  little 
difference  between  the  elastic  and  viscoelastic  response  for  small  time.  A 
comparison  of  an  undamaged  elastic  and  a  viscoelastic  response  Is  shown  In 
Figure  10  and  demonstrates  this  point.  An  Iterative  solution  Is  required  as 
discussed  earlier  since  damage  Is  occurring.  All  other  responses  are  simu¬ 
lated  for  a  viscoelastic  material  with  accumulated  damage.  The  damage  Is 
allowed  to  accumulate  as  long  as  the  circumferential  strain  Is  Increasing. 

When  the  circumferential  strain  decreases,  the  damage,  and  hence  E^,  Is  held 
constant. 

In  order  to  Implement  this  constant  damage,  two  loads,  P3  and  P^,  are 
added  to  load  P2.  The  response  to  P2  Is  viscoelastic  with  damage.  The  Ini¬ 
tial  conditions  are  found  from  the  state  that  exists  at  t^  due  to  P]^.  At  time 
t2»  the  negative  of  P2,  l.e.,  P3,  Is  applied.  Assuming  a  non-aging  material, 
the  Boltzman  superposition  principle  gives  the  response  to  P3  as  the  negative 
of  the  response  to  P2  delayed  by  x  -  t2  -  x^.  This  In  effect  "subtracts"  out 
the  damage  for  time  greater  than  t2«  If  a  step  decrease  In  load  to  zero  at 
time  t2  Is  desired,  no  further  loads  need  be  applied.  However,  the  desired 
loading  Is  a  linearly  decreasing  pressure  from  time  t2  to  13.  To  accomplish 
this,  the  load  P^,  which  Is  equal  to  P2,  Is  applied  at  time  t2.  The  response 
to  P4  Is  viscoelastic  with  constant  damage,  l.e.,  constant  E<]<.  This  Is  done 
so  that  we  may  "ramp  down"  the  load. 

The  linearly  decreasing  phase  of  the  loading,  l.e.,  t2  to  t3.  Is 
accomplished  by  applying  load  P5.  The  response  Is  viscoelastic  with  constant 
damage.  Load  Pg  Is  applied  at  time  t3  to  maintain  zero  load  for  time  greater 
than  t3.  The  response  to  Pg  Is  the  negative  of  the  response  to  P5  delayed  to 
X  ■  t3  -  t2. 


The  total  solution  Is  constructed  from  the  superposition  of  the  respon¬ 
ses  discussed  above  as  Indicated  below: 

1.  For  0  <  t  <  tj 

An  elastic  solution  with  damage.  An  Iterative  procedure  Is 
required. 


2.  For  ti  t  £  t2 

Viscoelastic  analysis  with  damage  for  the  step  load  P2.  Initial 
conditions  are  determined  from  the  state  at  t  ■  t^  from  1.  Solution 
form  Is  given  by  Eq.  (55).  An  Iterative  procedure  Is  required. 


I  ^  •  M  n 
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a.  The  solution  from  2.,  plus 

b.  The  negative  of  the  solution  from  2  shifted  by 
T  ■  t2  “  ti  hours,  plus 

c.  A  viscoelastic  solution  with  constant  damage  for  the  step  load 
P4.  Initial  conditions  are  determined  from  the  state  at 

t  >  t2  from  2.  Solution  form  Is  given  by  Eq.  (55),  plus 

d.  A  viscoelastic  solution  with  constant  damage  for  the  linearly 
varying  P5.  Initial  conditions  are  determined  from  the  state  at 
t  ■  t2«  Solution  form  Is  given  by  Eq.  (59). 

For  t3  ^  t 

a.  Solution  from  3.,  plus 

b.  The  negative  of  the  solution  from  3.c.  shifted  by 
T  ■  tq  -  t2  hours. 


VIII.  EXAMPLE  APPLICATION 


A  numerical  experiment  was  conducted  to  demonstrate  the  model  developed. 
The  solution  Is  for  the  time-dependent  circumferential  strain  of  a  filament- 
wound  composite  cylinder.  The  Internal  pressure  loading  Is  given  In  Figure  2 
with  Pq  -  5.516  MPq,  tj  -  0.00833  hr,  t2  -  0.025  hr,  and  t3  -  0.03333  hr. 

The  material  system  Is  a  Kevlar  49/HBRF-55A  filamentary  composite.  The 
Kevlar  49  fibers  are  considered  to  be  Isotropic  and  linearly  elastic.  Fiber 
properties  were  obtained  from  the  literature  [28].  The  reported  Young's 
modulus  and  Poisson's  ratio  are,  respectively,  131  GP^  and  0.20.  The  matrix 
material,  HBRF-55A,  Is  considered  to  be  Isotropic  and  linearly  viscoelastic. 

The  time-dependent  matrix  properties  used  here  were  determined  by  Hackett  and 
Dozier  [28]  from  creep  test  data.  The  expression  for  the  relaxation  modulus 
and  time-dependent  Poisson's  ratio  are  given  as,  respectively: 

Ej,(t)  -  152.5  exp  (-106. 5t)  +  40.23  exp  (-10.18t) 

+  154.6  exp  (-1.070t)  -  211.3  exp  (-0.092t) 

+  2346  exp  (-0.0075t) 

and 

Vt)  -  0.5  -  E„(t)/16548 
where  the  units  of  Ejj(t)  are  MP^. 

The  geometry  for  this  example  la  that  of  a  cylindrical  shell  having  an 
Inner  radius  of  38.10  mm  and  an  outer  radius  of  39.62  ran.  A  fiber  volume 
fraction  of  0.65  and  a  composite  failure  strain  of  0.0192  are  used.  The  fiber 
wrap  angle  (y)  is  +  55”.  Positive  wrap  angle  Is  shown  In  Figure  1. 

A  plot  of  the  simulated  response  with  and  without  damage  for  time  less 
than  t]^  Is  shown  In  Figure  11.  It  Is  seen  that  the  Inclusion  of  damage  In¬ 
duced  an  Increase  of  approximately  24%  over  the  undamaged  response.  This  Is 
Indicative  of  the  reduction  In  transverse  modulus  caused  by  the  damage.  The 
strain  reaches  a  given  level  In  smaller  time  for  the  damaged  composite. 

Implying  an  earlier  failure  than  without  damage. 

Figure  12  Is  a  comparison  of  the  modeled  response  with  and  without 
damage.  As  discussed  In  Chapter  VI,  the  strain  should  be  totally  recovered 
since  only  elastic  elements  In  the  model  are  damaged  and  no  plastic  behavior 
Is  Included.  Close  examination  of  Figure  12  shows  that  the  response  with 
damage  Indicates  a  very  small  negative  residual  strain.  This  Is  believed  to 
be  the  result  of  round-off  error  and  of  a  strong  sensitivity  to  Initial  con¬ 
ditions  that  was  observed  during  model  development.  Very  little  vlscoelastlc- 
Induced  strain  occurs  prior  to  time  t^.  The  damage  Is  mainly  due  to  the 
elastic  Increase  In  strain.  For  time  t]^  to  t2,  a  creep  under  the  constant 
load  Is  seen.  The  creep  predicted  for  the  damaged  case  Is  2.54  times  that 
predicted  for  the  undamaged  case.  The  damage  at  times  t]^  and  t2  are, 
respectively,  0:^^^  ■  0.7021  and  “  0.7386.  Thus,  the  creep  has  Induced  an 

additional  Increase  In  the  level  of  damage  of  approximately  3.65%. 
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Figure  11.  Comparison  of  damaged  and  undamaged  response 


for  t  < 


After  the  load  has  been  removed,  the  strain  Is  seen  to  recover  complete¬ 
ly  by  approximately  0.20  hours.  However,  the  damage  Is  not  recovered  and  the 
damage  present  at  time  t2  remains.  Using  Eq.  (69),  It  Is  calculated  that  the 
transverse  modulus  at  time  t2  retains  only  6.83Z  of  Its  Initial  value.  The 
majority  of  the  damage  Is  seen  to  occur  prior  to  time  tj^,  l.e.,  during  the 
linearly  Increasing  load,  with  a  smaller  portion  of  the  total  damage  occurring 
during  the  constant  load  phase.  This  Is  consistent  with  data  published  by 
Hahn  and  Kim  [38]  on  their  study  of  proof  testing  of  composite  materials. 

The  elastic  failure  load  (based  on  Instantaneous  values  of  the  moduli) 
are  calculated  at  Initial  time  and  for  maximum  damage.  Their  respective 
values  are:  “  9.977  and  ■  7.771  MP^^.  This  gives  a  22. IZ 

decrease  In  failure  load  as  a  result  of  tne  accumulated  damage.  The  calcu¬ 
lated  response  of  this  or  any  other  damaged  structure  to  future  loading  must 
be  based  on  the  effects  of  the  proof-test-induced  damage.  This  points  out  the 
Importance  of  considering  the  proof-test-induced  damage  In  determining  In- 
service  load  limits. 

Figure  13  contains  a  plot  of  the  circumferential  strain  versus 
pressure.  The  pressure  Is  plotted  as  the  vertical  so  that  a  qualitative  com¬ 
parison  between  Figure  13  and  Figure  7  can  be  made.  Comparison  Indicates  good 
agreement  and  gives  confidence  In  the  model.  The  creep  at  constant  load  and 
recovery  are  noted  by  the  change  In  strain  for  constant  pressure.  The  change 
In  stiffness  that  occurs  as  a  result  of  the  accumulated  damage  Is  Indicated  by 
the  change  In  slope  of  the  two  lines  drawn  tangent  to  the  response  curve.  The 
upper  line  corresponds  to  the  Initial  or  undamaged  stiffness  with  the  lower 
line  corresponding  to  stiffness  after  damage  has  accumulated. 


IX.  CONCLUSION 


A  Bodel  to  predict  the  vlscoelastlc/damage  response  of  a  filament-wound 
composite  cylindrical  pressure  vessel  to  a  proof-test  loading  was  developed. 
The  matrix  material  of  the  composite  system  was  assumed  to  be  Isotropic  and 
linearly  viscoelastic.  The  Kevlar  49  fiber  was  assumed  to  be  Isotropic  and 
linearly  elastic.  A  damage  model  which  relates  the  circumferential  strain  to 
a  damage  Index  was  developed.  This  damage  Index  was  In  turn  related  to  the 
transverse  modulus  as  a  first  approximation.  A  quadratic  relationship  between 
the  transverse  modulus  and  the  circumferential  strain  resulted.  An  Iterative 
solution  based  on  the  damaged  elastic  response  was  employed  to  predict  the 
time-dependent  circumferential  strain  during  the  linearly  Increasing  phase  of 
the  loading.  Thereafter,  an  Iterative  solution  based  on  the  elastlc-vlsco- 
elastlc  correspondence  principle  was  employed  to  predict  the  time-dependent 
circumferential  strain  %fhlch  Is  then  used  to  determine  the  transverse  modulus. 
A  new  expression  for  the  time-dependent  circumferential  strain  Is  then  found 
using  the  previously  calculated  transverse  modulus.  This  Iterative  procedure 
Is  continued  until  convergence  Is  achieved.  The  associated  elastic  solution 
resulting  from  use  of  the  elastlc-vlscoelastlc  correspondence  principle  Is 
Inverted  by  using  the  method  of  collocation  for  each  Iteration. 

The  simulation  model  developed  was  demonstrated  by  a  numerical  example. 
There  Is  good  qualitative  agreement  between  the  predicted  and  anticipated 
behavior.  After  unloading,  the  strain  was  totally  recovered.  However,  the 
transverse  modulus  has  retained  only  6.83X  of  Its  Initial  value.  The  response 
of  the  cylinder  to  future  loadings,  l.e..  In-service  loading,  must  Include  the 
effects  of  the  reduced  transverse  modulus. 

The  simulation  Is  well  founded  In  classical  viscoelastic  procedures  and 
the  damage  model  produces  results  which  are  consistent  with  reported  results 
and  should  prove  useful  In  future  work.  Future  efforts  should  Include  experi¬ 
ments  to  verify  the  model  developed,  as  well  as  efforts  to  Investigate  the 
effect  on  the  response  of  different  straln.-damage-transverse  modulus  relation- 
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APPENDIX.  COMPUTER  LISTING  OF  SOLUTION  MODEL 


PROGRAM  TO  CALCULATE  THE  VISCOELASTIC  HOOP  STRAIN 
OF  A  FILAMENT  WOUND  CYLINDRICAL  PRESSURE  VESSEL 
INCLUDING  A  “DAMAGE"  EFFECT. 


BY  JOHNNY  L.  PRATER  UNDER  DIRECTION  OF  DR.  ROBERT  M.  HACKETT  FOR 
PARTIAL  FULFILLMENT  OF  REOUIRMENTS  FDR  PH.D.  DEGREE  IN  SOLID 
MECHANICS  AT  THE  UNIVERSIT  OF  ALABAMA' IN  HUNTSVILLE.  1984-1985. 

LOAD-  M*T 

PARAMETERS 

EF  -  YOUNG'S  modulus  OF  THE  FIBER  (PSD 
PF  -  POISSON'S  RATIO  OF  THE  FIBER 
VF  -  FIBER  VOLUME  FRACTION 
VM  -  MATRIX  VOLUME  FRACTION 

EMI  -  INITIAL  RELAXATION  MODULUS  OF  MATRIX  (PSD 
PMI  -  INITIAL  POISSON'S  RATIO  OF  MATRIX 

GAMMA  -  ANGLE  BETWEEN  THE  FIBER  LONGITUDINAL  DIRECTION  AND  THE 
CYLINDER  THETA(  HOOP  )  DIRECTION  (DEG). 

RIN  -  INSIDE  RADIUS  OF  CYLINDER  (INCH) 

ROUT  -  OUTER  RADIUS  OF  CYLINDER  (INCH) 

XLOAD  -INTERNAL  PRESSURE  (PSD  OF  THE  CYLINDER 

THK  -  COMPOSITE  FAILURE  STRAIN 

FLI  -  ELASTIC  FAILURE  LOAD  (PSD 

EPTHTH  -  HOOP  STRAIN 

EPMM  -  MERIDIONAL  STRAIN 

EPRR  -  RADIAL  STRAIN 

PSTRN  =  HOOP  STRAIN  IN  LAPLACE  -DOMAIN 
PLOAD  =  LAPLACE  TRANSFORM  OF  THE  LOAD 
T  -  TIME  (HOURS) 

NLOOPl  »  NUMBER  OF  LOOPS  USED  FOR  CONVERGENCE  OF  ELASTIC 
PHASE  (  O  <  T  (  T-1  ), 

NL0QP2  -  NUMBER  OF  LOOPS  USED  FOR  CONVERGENCE  OF  VISCOELASTIC 
PRASE  (  T1  (  T  ) . 


IMPLICIT  REAL *8  <A-H.P-Z) 

REAL*8  LOAD,  LAMBDAI.  LAMBDA 

DIMENSION  STRN<200)  . PSTRN (6) . GV (6)  .  STN ( 6) .5(6) .V(6> ,A(&.6) 
DIMENSION  GVl (6) .QV2(6) .AA(6.6) .STNl (6) .5TN2(6) 

DIMENSION  EMP(6) .PMP(6) .GVCD(6> .STNCD(6) 

0PEN(UNIT-7, FILE-' VISD. OUT' . STATUS- ' NEW ' > 
OPEN(UNIT-14.FILE-'VISD.PLT' .STATUS-’NEW’ ) 

0PEN(UNIT-15. FILE-' VISD. PRS' . STATUS-' NEW' ) 


WRITE(6.*) 

WRITE(6.*)'  THE  SOLUTION  WILL  BE  FOR  AN  ELASTIC  COMPOSITE' 
WRITE(6.*)'  CYLINDER  WITH  A  LOAD  AS  SHOWN  BELOW.' 
WR1TE(6.*) 


A-1 


nnn 


(Tl. XLDAD) • 


WRITE(6. •) ’ 

WRITE(6.*>’  I  *.*•***•  <T2.XL0AD)‘ 

WRITE(6.*>'  I  •  •  • 

WRITE<6.'*)*  I  •  •! 

WRITE<6,*)*'  I  •  •• 

WRITE(6.*)’  I*  •• 

WRITE<6,*)*  • _ • _ • 

MRITECb.*)*  (0.0)  (T3.0>* 

WRITE (6. *> 

WRITE (6. •) 

WRITE(6. *> 

XLOADaSOO. 

T 1 - . 008333333333 
T2=.025 

T3-. 033333333333 
GAMMA=5S. 

RIN>1.S 

ROUT-I.SA 

WRITE<7.*)'  LOAD  (PSD  -  • . XLQAD 

WRITE(7,*>’  Tl  (HR)  -  » .Tl 

WRITE<7.*)*  T2  (HR)  -  •.T2 

WR1TE(7.*)’  T3  (HR)  -  •.T3 

WRITE{7.*)'  GAMMA  (DEG)  =  ’.GAMMA 

WRITE(7.*)*  RIN  (IN)  -  • .RIN 

WRITE(7.*)*  ROUT  (IN)  »  ’.ROUT 

WRITE (6.*)  ’  NUMBER  OF  LOOPS  FOR  ELASTIC  SOLUTION  ??7 
REAO(S. *)NL0QP1 

WRITE (6.*)’  NUMBER  OF  LOOPS  FOR  CREEP  PHASE  ???  ’ 

READ (3. •)NL00P2 

w»o. 

N0«O 

SLOPEl-XLOAD/Tl 


C  FAILURE  STRAIN 

THK-.0192 


TIME-DEPENDENT  MATRIX  PROPERTIES  WERE  DETERMINED  BY  LEAST  SQUARES 
CURVE  FIT  TO  CREEP  TEST  DATA  BY  JAN  DOZIER. 

INITIAL  VALUE  OF  RELAXATION  MODULUS  AND  POISSON’S  RATIO  OF  HPHF-30A 

EMI -359933. 

PMI-  0.3501 

WRITE(7.*)’  EMI  -  ’ .EMI 
WRITE(7.*)’  PMI  -  ’.PMI 

C  YOUNG’S  MODULUS  AND  POISSON’S  RATIO  OF  XEVUAR  49  FIBER. 

EF-  19000000. 

PF-  0.20 

WRITE (7.*)’  EF  -  ’.EF 
WRITE(7, •) ’  PF  -  ’ .PF 


C 


FIBER  VOLUME  FRACTION 


onn  nn  nnn 


VF=  0.65 

WRITE <7. •> •  VF  =  V.  VF 


ELASTIC  RESPONSE  WITH  DAMAGE 
O  <  T  <  T1 


PI-1. 

PI-  A. *ATAN(PI) 

NTIMES-100 
DELT-Tl/NTIMES 
GAMMA-GAMMA •P I / 1 80. 

DO  2000  IX-l.NLOQPl 
T-DELT 

DO  1000  I-l.NTIMES 
LOAD-SLOPE 1*T 

CALCULATE  TRANSVERSE  MODULUS  RESULTING  FROM  DAMAGE  CALCULATED 
IN  LAST  LOOP. 

IF<rX  .NE.  1  )ET-ETI*<STRN(I)/THK-1)**2 

CALL  SOLN ( GAMMA . R I N . ROUT . VF . EF . PF , EM I . PM I , THK , LOAD , 

•  FLI.EPTHTH.EPMM.EPRR.U.ND.ET.W) 

STRN  < I ) -EPTHTH 
TL0G«L0G10<T> 

IF < I X . EQ. NLOOPl > THEN 

PRESS-SL0PE1*T*689S. / 1000000. 

WRITE <14, 101)T,STRN<I) 

WRITE!  IS. 101) STRN <l> .PRESS 

END  IF 

101  F0RMAT<2aJC.E14.S)  X 

T-T+DELT 
1000  CONTINUE 

IFdX.EQ.  DTHEN 
ETI=ET 

WRITE(6.*)’  FAILURE  LOAD  -  '•.FLI 
WRITE!?,*)’  FAILURE  LOAD  -  ’.FLI 

END  IF 
ND-1 

WRITE!6.*)’  STRN  -  ’,STRN!100> 

2000  CONTINUE 

STRN I -STRN  <  NT IMES ) 

W1 -STRN I /THK 


CONSTANT  LOAD  PHASE 
! CREEP  WITH  DAMAGE) 


CALCULATE  THE  LAPLACE  PARAMETERS. 


DO  201  I>1.6 

S<I)*«EXP(RR*I-Q> 

V(I>«1/S<I) 

201  CONTINUE 

ETP=ETI* (Wl-1) **2 
DO  2100  IX-1.NL00P2 
DO  1100  I«1.6 


C  LAPLACE  TRANSFORM  OF  XLOAD 

PLOAD=XLQAD/S<I) 

C  LAPLACE  TRANSFORM  OF  MATRIX  PROPERTIES. 

EMP<I)-S<I)*(22120. /<S<I>+106.3>+5835./<S<I)'HO. 10) 

•  +22420. /<S(I)-H.07>-30&40./<S(I)*. 092) 

•  +340200./ (S(I) +.0075) > 

PMP  < I ) ». 5-EMP  < I ) / 2400000. 

C********************************************.*****..****..*.*; 

C  CALCULATE  LAPLACE  TRANSFORM  OF  DAMAGED  TRANSVERSE  MODULUS 

IFdX  .QT.  DTHEN 

ETB- < 1 /S  < I ) -2*PSTRN< I ) /THK  +  ( AC*  *2) /S ( I) /THK • *2) 

DO  310  IJ-1.& 

ETB=ETB+2*AC*QV<I J)/ <S(I)+1/V< IJ) > /THK**2 
DO  320  IK-1.6 

ETB-ETB+GV<I J) *60 < IK4 / <S < I ) +1 /V ( I J> +1/V ( IK > ) /THK**2 
320  CONTINUE 

310  CONTINUE 

ETP-ETB*S<I) -ETI 

ENDIF 


CALL  SOLN<QAMMA.RIN,ROUT.VF.EF.PF.EMP(I> .PMPU) .THK.PLOAD. 

•  PFL.PSTRN(I) .PEPMM.PEPRR.PU.ND.ETP.Wl) 

C  ••••  RIGHT  HAND  SIDE  VECTOR  TO  USE  IN  SOLN  FOR  COEFFICIENTS  <GV»S) 

STN<I)-(PSTRN<I)-  STRNI/Sd)  )*Sd) 

1100  CONTINUE 
ND-1 


C*************************************^**. . . . 

C  PERFORM  LAPLACE  TRANSFORM  INVERSION  AND  SOLVE  FOR  CONSTANTS 

C  TO  BE  USED  IN  TIME  DEPENDENT  MAXIMUM  HOOP  STRAIN. 

DO  220  1-1,6 
DO  220  J-1.6 

Ad.  J)-<Vd)*V<J>/(Vd)+V<J)  )-Vd>>/Vd  ) 

220  CONTINUE 


C 


•••••  assemble  MATRIX  FO  INVERSION 


n  n 


CALL  SOLVe(A.STN.&) 


AC=STRNI 
DO  240  I>>1.6 
QV(I)»STN<I) 
AC»AC-GV{1) 
240  CONTINUE 


C  CALCULATE  TIME  TO  FAILURE  <USED  TO  CHECK  CONVERGENCE) 


TNI -50. 

DO  1=1,10 

F1=AC-THK 

DF1=0. 

DO  32  IJ=1.6 

Fl«Fl-*-GV<I  J)  •EXP<-TN1/V<I  J)  ) 
DF1-DF1-GV<1J>*EXP<-TN1/V(1 J>  >/V<l J> 

32  CONTINUE 

TN2-TN1-F1/DF1 
TERR-TN2-TN1 
TN1-TN2 
31  CONTINUE 

WRITE<6.*)’  T  FAIL  -  ’.TN1.»  TERR  -  ‘.TERR 
2100  CONTINUE 


C  CALCULATE  STRAIN  AT  TIME  T2  <USED  AS  INITIAL  CONDITION  FOR 

■■^RNIl-AC 
T»  •2-Tl 
DO  312  1-1.6 

QTRNII-STRNI1*QVU>  -EXPl-T-SlI)  ) 

312  CONTINUE 


C 


ASSOCIATED  SOLUTION 
FOR  RAMP  DOWN  PHASE 


C  CALCULATE  INITIAL  PARAMETERS  FOR  RAMP  DOWN 

ND-0 

W2-STRN11/THK 
SLOPE— XLOAO/  1T3-T2) 

ET=ETI-<W2-l>-*2 

CALL  SOLN <GAMMA. RIN. ROUT, VF.EF.PF. EMI, PMl .THK. XLOAD. 
•  FH .EPTHTH.EFMM.EPRR,U.ND.ET.W2) 

SRATE— EPTHTH/  <  T3-T2) 

C  "ASSOCIATED  ELASTIC"  SOLUTION 


A- 5 


DO  1200  I-I.6 


c 


LAPLACE  TRANSFORM  OF  THE  LOAD 


PLOAD=SLOPE/S ( I) **2 

CALL  SOLN  IGAMMA .  R IN .  ROUT .  VF .  EF .  PF‘.  EMP  (  I  >  .  PMP  ( I  >  .  THK  ,  PLOAD . 
•  PFL . PSTN . PEPMM .PEPRR.PU.ND.ET.W2) 

C  RIGHT-HAND  SIDE  VECTOR  TO  USE  IN  SOLN  FOR  COEFFICIENTS 

STNl (I )=(PSTN-  SRATE* (V ( I ) **2) ) *S( I > 

STN2  <  I  >  =-STN.l  <  I ) 

1200  CONTINUE 


C  DEVELOPE  MATRIX  OF  LAPLACE  PRAMETERS  (  A  > 

CALL  AMAT(V.0.0.A.6) 

C  SOLVE  FOR  COEFFICIENTS  OF  EXPONENTIAL  TERMS. 

DO  11  1-1.6 
DO  11  J-1.6 
AA(I. J)-A(I. J) 

11  CONTINUE 

CALL  SOL VE< A. STNl, 6) 

CALL  SOLVE<AA.STN2.6> 

ACl-O. 

AC2-0. 

ADI -SPATE 
AD2— SRATE 
DO  2S0  1-1.6 
ACl-ACl-STNl <I) 

AC2-AC2-STN2  < 1 ) 

ADl-ADl-^SXNl  (I)  /VTI) 

AD2-AD2+STN2  < I ) /V ( I ) 

GVl <I)=STN1 (I) 

QV2<I)=STN2<I) 

250  CONTINUE 


C  *•••  CONSTANT  DAMAGE  CREEP  ••• 


ND-1 

ETP-ETI*  <W2-1) ••2 

CALL  SOLNCGAMMA. RIN. ROUT. VF.EF.I^F. EMI. PMI .THK , XLOAD. 

•  FLI.EPTHTH,EPMM.EPRR.U.ND.ETP.W2> 

STRNI2-EPTHTH 

DO  1300  1-1,6 
PLOAO-XLOAD/S(I) 

CALL  SOLN ( GAMMA , R I N . ROUT; . VF . EP . PF . EMP ( I > . PMP (I). THK. PLOAD. 

•  PFL. PSTRN<1> .PEPMM, PePRR,PU.ND,ETP,W2) 

C  ••••  RIGHT  HAND  SIDE  VECTOR  TO  USE  IN  SOLN  FOR  COEFFICIENTS 
8TNCD(  I>»(PSTRN(I)  -  STRN12/S«I)  )  -Sd) 


A- 6 


PERFORM  LAPITACE  TRANSFORM  INVERSION  AND  SOLVE  FOR  CONSTANTS 
TO  BE  USED  IN  TIME  DEPENDENT  MAXIMUM  HOOP  STRAIN. 

DO  222  I>1.6 
DO  222  J=1.6 

Ad.  J)  =  <V(I)  •VC  J)  /<V<  I)  ♦V<  J>  >~V<I>  >  /V<I> 

222  CONTINUE 

••••  ASSEMBLE  MATRIX  FOR  INVERSION. 

CALL  S0LVE(A.STNC0.6> 


ACCD=STRNI2 
DO  244  1-1.6 
GVCD( I) -STNCDCI) 
ACCO>ACCO-QVCD< I ) 
244  CONTINUE 


CALCULATE  CIRCUMFERENTIAL  STRAIN  IN  TIME  DOMAIN. 

T-.OOl 

DO  300  LL-1.12S 
STRNT-0. 

DO  3S0  1-1,6 

STRNT-STRNT+GVdJ  •EXP<-T»Sd)  > 

I F  < T. QT . T2-T I ) STRNT-STRNT-QV < I ) -EXP < - (T- ( T2-T I > ) •S ( I >  > 

I F  <  T . QT . T2-T 1 >  STRNT-STRNT+QVCD ( I > •EXP  < - <  T- ( T2-T 1 >  > •S ( I > ) 

I F  < T. QT . T2-T I ) STRNT-STRNT+QVl ( I ) •EXP <- ( T- {T2-T I >  > •S ( I >  > 

IF<  T  .QT.  T3-T1)STRNT-STRNT+QV2(I>^EXP<-(T-<T3-T1) )^S(I> ) 
350  CONTINUE 

STRNT-STRNT+AC 
I F  <  T . GT . T2-T I ) STRNT-STRNT-AC 
I F  <  T . GT . T2-T I ) STRNT-STRNT+ACCD 
I F  <  T . GT . T2-T I ) STRNT-STRNT*AC1 ♦AD  <  T- ( T2-T I >  > 

IF<T  .QT.  T3-TI  )  STRNT-STRNT+AC2  ♦AD2* <T- < T3-TI ) ) 

PRESS=XL0AD«6e95. /lOOOOOO. 

IF  CT, QT. T2-TI ) PRESS-PRESS^ (SLOPED (T- tT2-Tl ) > ) •6895. / lOOOOOO. 
IF  <  T . GT . T3-T 1 ) PRESS-0. 

TT-T+Tl 

TLOG-LOQIO(TT) 

WRITE  (14.  lODTT.STRNT 
WRITE  (15.  lODSTRNT. PRESS 
T-T+.0015 
300  CONTINUE 


SUBROUT I NE  SOLN  <  GAMMA , R I N . ROUT . VF . EF . PF . EM I . PM I . THK . LOAD . 

FL 1 . EPTHTH . EPMM . EPRR .U.ND.ETl.W) 

IMPLICIT  REAL*8  <A-H.P-Z> 

REAL*8  LAMBDA. LOAD. LAMBDA I 
PI-I. 

PI=4. *ATAN(PI) 

S1=SIN<GAMMA> 

Cl=COS< GAMMA) 

S2=S1**2 

C2=C1**2 

S4=S1**4 

C4-C1*«4 

VM-  1.0  -  VF 

GF=  EF/ <2. • <1.+PF) ) 

GMI=  EMI/ (2. • (l.+PMI) ) 

ELI-  VF*EF+VM*EMI 

IF  <ND  .EQ.  DGO  TO  2 

ETI - < 1 . -W>  **2/ ( VF/EF+VM/EMI ) 

PLTI-  VF*PF-*-VM*PMI 
PTLI-  PLTI-ETI/ELI 

PTTI-  VF*PFWM*PMI*  (  <1+PMI-PLTI*  (EMI/ELI)  >/ 
<l-PMI**2+PMl*PLTI«<EMI/ELt) ) ) 

BNI-  (GF*  <PI+4*VF)-*-GMl*  <PI-4*VF)  >/ 

■  <QF*  (PI-4*VF)4-GMI*  (P1+4*VF)  ) 

QLTI-  (GMI/2) *< ( <4-PI)+PI*BNI)/4+4*BNI/ ( (4-Pn •BNI+PI )  ) 

LAMBDAI -  1 / <  < 1 ♦PTT I ) • ( 1 -PTT I -2*PLT 1 -PTL I >  > 

CLLl-  <1-PTTI*«2) *LAMBDAI*ELI 

CTTI-  < I -PLTI ♦PTLI) -LAMBDAI •ETI 

CLTI-  PTLI-<l+PTTI) -LAMBDAI-ELl 

CTR I -  (PTTI -PLT I -PTL I ) -LAMBDAI -ET I 

GTR I -  < I -PTT I -2-PLT I -PTL I ) -LAMBDA I -ET I /2 

CRRI-  CTTI 

CMMI-S4-CLLI+C4-CTTI+2-S2-C2- (CLTI-2-GLTI ) 
CTHTHI-C4-CLLI+S4-CTTI+2-C2-S2-<CLTI-2-GLTI) 

CMTH I - <  S4+C4 ) -CLT I* ( CLL I -CTT I -4 -GLT I ) •S2-C2 

CMRI-S2-CLTI+C2-CTRI 

CTHRI -C2-CLT I +S2-CTR I 

ALPHA2-CTHTH I /CRR I 
ALPHA-SORT  <ALPHA2) 

BET A- < CMTH I -CMR I ) / CRR I / ( 1 -ALPHA2 ) 


CALCULATE  CAPA,  CAPB.  AND  EP»W  ,  NEEDED  TO  CALCULATE  U  (DISPLACEMENT) 

All-2- (ROUT-- < I -ALPHA) -RIN-- ( 1+ALPHA) ) • (CMTHI -ALPHA-CMR I ) 

/ (1+ALPHA) /R1N--2 

A12-2- (ROUT-- (I -ALPHA) -RIN--< 1-ALPHA) 5 • (CMTHI-ALRHA-CMRI ) 

/ <1-ALPHA) /RIN--2 

A13-<ROUT--2-RIN--2) • (CMMl-BETA- (CMTHI +CMRI ) ) /RIN--2 
A21- (CTHRI -ALPHA -CRR I ) /RIN-- ( 1-ALPHA) 

A22-(CTHRI-ALPHA-CRRI ) /RIN-- ( l-ALPHA) 

A23-CMR I -BETA- ( CTHR I -CRR I ) 

A31 - ( CTHR I -ALPHA-CRR I ) /ROUT— ( 1 -ALPHA) 

A32-(CTHR1-ALPHA-CRR1 ) /ROUT-- ( 1 -ALPHA) 

A33-CMR I -BETA • ( CTHR 1 -CRR I ) 


B1 1-A22*A33-A32*A23 
B12—  <A21  •A33-A31  •A23) 

B13>A21 *A32-A22*A31 
B21  — (A12*A33-A32*A13) 

B22=A1 1 •A33-A31 *AJ3 
B23=- (A1 1 *A32-A31 •A12) 

B31=A12*A23-A22*A13 
B32—  (All  •  A23-A2 1  •  A 1 3 ) 

B33=A1 1 •A22-A12*A21 

DELTA=A1 1*A22*A33+A12*A23*A31+A13*A21*A32 
•  -A13*A22*A31-A12*A21*A33-AH*A32*A23 


CAPA^LOAO* (B1 1-B21 ) /DELTA 
CAPB=LQAD*  <B12-B22) /DELTA 
EPMM-LOAD*  <B13-B23) /DELTA 


CALCULATE  RADIAL  DISPLACEMENT  (  U  >.  HOOP  STRAIN  (  EPTHTH  > 
AND  RADIAL  STRAIN  (  EPRR  ). 


U»CAPA*R**ALPHA+CAPB/R**ALPHA+BETA*EPMM*R 


CALCULATE  THE  FAILURE  LOAD 

FLI=LOAD*THK*R/U 

EPTHTH-U/R 

EPRR-CAP A • ALPHA • R • * ( ALPHA- 1 ) -CAPB* ALPHA / R • • <  ALPHA* 1 > 
•  +BETA*EPMM 

RETURN 
END 


SUBROUTINE  SOLVE (E.G.N) 

IMPLICIT  REAL*a  <A-H.P-Z> 

DIMENSION  E<N.N).Q(N) 

C  SIMPLE  GAUSS  ELIMINATION 

C  FORWARD  ELIMINATION 

DO  200  J-l.N-1 

IF<E( J. J) .NE.O. »GO  TO  205 

write (&•*)’ ******  E(J.J)  IS  ZERO  FOR  J  -  '•J 
RETURN 

205  EJJ-  E(J.J) 

DO  201  I>J.N 
E(JrI)-E<J.I)/EJJ 
201  CONTINUE 

G< J>=Q< J) /EJJ 
DO  220  JJ=J*1.N 
CJ=E<JJ. J) 

DO  230  I=J.N 

E<JJ. I)-E(JJ.I>-E(J.I)*CJ 
230  CONTINUE 

Q<JJ)-Q<3J)-Q<J)*CJ 
220  CONTINUE 
200  CONTINUE 


C  BACK  SUBSTITUTION 


f 


Q(N>«Q<N)/E(N.N> 

DO  240  J-N-l.l.-l 
DO  250  K-J+l.N 
Q< J)-Q< J>-G(K>*E< J.K) 
250  CONTINUE 
240  CONTINUE 
RETURN 
END 


SUBROUTINE  AMAT (V. TD. A. N) 

REAL*a  V.TD.A 
DIMENSION  V(N).A(N.N> 

DO  220  IW-1.6 
DO  220  IV«1.6 

A(IW.  IV)s(V(  IV)  •0<IW)  •EXP<TO/V(IV)  )  /  ( V  ( IV> '^V  ( IM>  > 

•  ♦V<IW)»*2/V<IV)-V(IW) 

•  -V(IW)«TD/V<IV>)/V(IW) 

220  CQNT I NUE 

RETURN 

END 
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